In [1]:
import sklearn
import numpy as np
from scipy import ndimage

from matplotlib import pyplot as plt
%matplotlib inline
import matplotlib.patches as mpatches

import skimage
from skimage import io
from skimage.filter import threshold_otsu

from skimage.feature import peak_local_max
from skimage.morphology import watershed
In [2]:
img1 = skimage.io.imread("img/cyano1.jpg")

bplt, (sub1) = plt.subplots(nrows=1,ncols=1,figsize = (9,9) )
sub1 = plt.imshow(img1)

arrow_args = dict(arrowstyle="->")

sub1.axes.annotate('Heterocyst'
                   ,xy= (1600,700)
                   ,xytext = (700,300)
                   ,xycoords = 'data'
                   ,textcoords = 'data' #'offset points'
                   ,arrowprops = arrow_args
                   ,fontsize = 14
                   )
sub1.axes.annotate('Vegetative'
                   ,xy = (1080,960)
                   ,xytext = (400,600)
                   ,xycoords = 'data'
                   ,textcoords = 'data' #'offset points'
                   ,arrowprops = arrow_args
                   ,fontsize = 14
                   )
Out[2]:
<matplotlib.text.Annotation at 0x1500b860>
In [3]:
def boundbox1(tryit):
#     np.nonzero(np.apply_over_axes( np.argmax, thelabel,[1]))[0].min()
    minX = np.nonzero(np.apply_over_axes(np.argmax,(tryit),[1]))[0].min()
    maxX = np.nonzero(np.apply_over_axes(np.argmax,(tryit),[1]))[0].max()
    minY = np.nonzero(np.apply_over_axes(np.argmax,(tryit),[0]))[1].min()
    maxY = np.nonzero(np.apply_over_axes(np.argmax,(tryit),[0]))[1].max()
    return [minX,maxX,minY,maxY]
In [4]:
def buildrect(arr):
    
    minY,maxY,minX,maxX = arr
    
    rect1 = mpatches.Rectangle(xy=(minX,minY)
                           ,width = maxX - minX
                           ,height = maxY - minY
                           ,fill = False
                           ,edgecolor = 'red'
                           ,linewidth = 2)
    return rect1
In [5]:
def img_to_labels1(inp_img):
    
    thresh = threshold_otsu(inp_img)
    otsu0 = inp_img < thresh
    dist0 = ndimage.distance_transform_edt(otsu0)
    j=34
    peaks0 = peak_local_max(dist0 ,indices = False
                            ,footprint=np.ones((j,j))) 
    print peaks0.cumsum()
    marker0 = ndimage.label(peaks0)[0]
    labels0 = watershed(-dist0, marker0, mask=otsu0)
    
#     plt.imshow(labels0)
#     plt.show()
    return labels0
In [6]:
def make_windows(inp_labels, inp_img):
    bb0 = []
    for i in range(1,inp_labels.max() + 1):
        label_i = (inp_labels == i)
        bb0.append(boundbox1(label_i))

    rect_arr = [buildrect(bb0[i]) for i in np.arange(bb0.__len__())]

    fig = None
    s0 = None
    fig, (s0) = plt.subplots(nrows = 1, ncols = 1, figsize = (12,12))
    s0.imshow(inp_img, cmap='gray')
    for rect in rect_arr:
        s0.add_patch(rect)
    fig.show()
    return fig
In [7]:
def algo1(inp_img):
    labels = img_to_labels1(inp_img)
    out_fig = make_windows(labels,inp_img)
    return out_fig
In [8]:
img1_gray = skimage.color.rgb2gray(img1)
In [10]:
plt.imshow(img1_gray, cmap='gray')
Out[10]:
<matplotlib.image.AxesImage at 0x18880588>
In [11]:
cross = img1[1200:2000,0:500]
In [12]:
cross = skimage.color.rgb2gray(cross)
In [14]:
algo1(cross)
[ 0  0  0 ..., 25 25 25]

Out[14]:
In [15]:
inp_img = cross
thresh = threshold_otsu(inp_img)
otsu0 = inp_img < thresh
dist0 = ndimage.distance_transform_edt(otsu0)
j=34
peaks0 = peak_local_max(dist0 ,indices = False
                        ,footprint=np.ones((j,j))) 
print peaks0.cumsum()
[ 0  0  0 ..., 25 25 25]

In [16]:
img = peaks0
peaks_array =[]

peak_counter = 0

for y in range(0,img.shape[0]):
    for x in range(0,img.shape[1]):
        
        if (img[y,x] == True):

            peak_counter = peak_counter + 1
            peaks_array.append([peak_counter,y,x])
            
In [17]:
img = cross
offset = 40

bplt, (sub1) = plt.subplots(nrows=1,ncols=1,figsize = (25,25) )
sub1 = plt.imshow(img, cmap='gray')

arrow_args = dict(arrowstyle="->")

for i in peaks_array:
    
    sub1.axes.annotate('Peak '+ str(i[0] - 1)
                   ,xy= (i[2],i[1])
                   ,xytext = (i[2]+offset,i[1]+offset)
                   ,xycoords = 'data'
                   ,textcoords = 'data' #'offset points'
                   ,arrowprops = arrow_args
                   ,fontsize = 12
                   )
In [18]:
X = np.array(peaks_array)[:,(1,2)]
In [19]:
from sklearn.neighbors import kneighbors_graph
In [20]:
weighted_graph = kneighbors_graph(X, 5, mode='distance')
In [21]:
def strand_v1(dist, pres, past):
    """
    args:
        dist: sklearn distance kneighbors_graph
        pres: present label
        past: past label
    """
    m, n = dist.shape
    n1nbr = None
    for i in range(n):
        if i in past or dist[pres, i]==0:
            continue
        if n1nbr is None or dist[pres,i]<=dist[pres,n1nbr]:
            n1nbr=i
#     print pres, past, n1nbr
    if len(past) == 0 or n1nbr is not None:   
        return strand_v1(dist, n1nbr, past+[pres])
        
    return past + [pres]
In [26]:
def strand_v2(dist, pres, past):
    """
    Uses a mean std deviation metric to ensure neighbors are cells reasonably close
    to one another, based on 
    args:
        dist: sklearn distance kneighbors_graph
        pres: present label
        past: past label
    """
    m, n = dist.shape
    n1nbr = None
    for i in range(n):
        if i in past or dist[pres, i]==0:
            continue
        if n1nbr is None or dist[pres,i]<=dist[pres,n1nbr]:
            n1nbr=i
#     print pres, past, n1nbr
    if len(past) == 0 or n1nbr is not None: 
        rawr = dist[]
        if n1nbr <  and np.std(past+[pres])==0:
            return strand_v2(dist, n1nbr, past+[pres])
        
    return past + [pres]
In [185]:
rawr = np.random.rand(10)*500

print rawr

for i in range(1, len(rawr)):
    print mean_std(rawr[:i])
[ 407.3219218   458.21728889  148.86451414   68.02138391  378.54867983
   91.91384786  110.33535318  478.40299549  252.90944502  476.76539288]
(407.32192180083103, 0.0)
(432.76960534446505, 25.447683543634042)
(338.13457494190351, 135.43744321178494)
(270.60627718347877, 165.64325752127829)
(292.19475771373106, 154.31915297923567)
(258.81460607068101, 159.42555242313685)
(237.60328422915831, 156.47714936402963)
(267.70324813617725, 166.63287724630794)
(266.05949223434266, 157.17176218128432)

In [29]:
strand_v2(weighted_graph, 3, [])
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-29-8e40da581a86> in <module>()
----> 1 strand_v2(weighted_graph, 3, [])

TypeError: strand_v2() takes exactly 4 arguments (3 given)
In [30]:
strand_v2(weighted_graph, 19, [])
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-30-517d77a771e4> in <module>()
----> 1 strand_v2(weighted_graph, 19, [])

TypeError: strand_v2() takes exactly 4 arguments (3 given)
In [188]:
strand_v1(weighted_graph, 15, [])
Out[188]:
[15, 16, 14, 11, 17, 18, 20, 21, 22, 23, 24, 19, 13, 12]
In [142]:
strand_v1(weighted_graph, 5, [])
Out[142]:
[5, 6, 4, 2, 0, 1, 3, 7, 9, 10, 12, 13, 19]

Test Polynomial Fitting

In [31]:
hand_labeled = [19, 13, 12, 10, 9, 7, 8]
# labels=tuple(strand_v1(weighted_graph, 19, []))
labels = hand_labeled
print labels
[19, 13, 12, 10, 9, 7, 8]

In [32]:
X[labels,:]
Out[32]:
array([[501, 480],
       [459, 442],
       [458, 442],
       [428, 393],
       [427, 392],
       [410, 344],
       [416, 287]])
In [33]:
coefficients = np.polyfit(X[labels,1], X[labels,0], len(labels)-1)
print coefficients
[  2.81586458e-11  -4.21901677e-08   1.79701007e-05   2.46024568e-03
  -4.01857258e+00   1.10485572e+03  -9.93364666e+04]

c:\Users\louiery\AppData\Local\Continuum\Anaconda\lib\site-packages\numpy\lib\polynomial.py:587: RankWarning: Polyfit may be poorly conditioned
  warnings.warn(msg, RankWarning)

In [34]:
polynomial = np.poly1d(coefficients)
In [35]:
ys = polynomial(X[labels,1])
In [36]:
fig,ax=plt.subplots(figsize=(15,15))
ax.imshow(img, cmap='gray')
ax.plot(X[labels,1], ys, 'r', lw=3)
ax.set_ybound(img.shape[0],0)
ax.set_xbound(img.shape[1],0)

Most Likely Path Finding with Polynomial Derviatives

In [37]:
img = cross
offset = 40

bplt, (sub1) = plt.subplots(nrows=1,ncols=1,figsize = (25,25) )
sub1 = plt.imshow(img, cmap='gray')

arrow_args = dict(arrowstyle="->")

for i in peaks_array:
    
    sub1.axes.annotate('Peak '+ str(i[0] - 1)
                   ,xy= (i[2],i[1])
                   ,xytext = (i[2]+offset,i[1]+offset)
                   ,xycoords = 'data'
                   ,textcoords = 'data' #'offset points'
                   ,arrowprops = arrow_args
                   ,fontsize = 12
                   )
In [40]:
# east_strand = strand_v1(weighted_graph, 19, [])
east_strand = hand_labeled
In [41]:
north_bound = east_strand + [11, 6]
print north_bound
[19, 13, 12, 10, 9, 7, 8, 11, 6]

In [42]:
west_bound = east_strand + [11, 14]
print west_bound
[19, 13, 12, 10, 9, 7, 8, 11, 14]

In [43]:
south_bound = east_strand + [11, 17]
print south_bound
[19, 13, 12, 10, 9, 7, 8, 11, 17]

In [44]:
def show_path_and_der(coords, labels, do_plot=True):
    labels = tuple(labels)
    coefficients = np.polyfit(coords[labels,1], coords[labels,0], len(labels)-1)
    polynomial = np.poly1d(coefficients)
    ys = polynomial(coords[labels,1])
    dcoefficients = np.polyder(coefficients)
    print dcoefficients
    dpolynomial = np.poly1d(dcoefficients)
    dys = dpolynomial(coords[labels,1])
#     d2coefficients =np.polyder(coefficients, 2)
#     print d2coefficients
#     d2polynomial = np.poly1d(d2coefficients)
#     d2ys = d2polynomial(coords[labels,1])
    
    if do_plot:
        fig,ax=plt.subplots(figsize=(10,10))
        ax.imshow(img, cmap='gray')
        ax.plot(coords[labels,1], ys, 'r', lw=3)
        ax.plot(coords[labels,1], dys, 'g', lw=3)
#         ax.plot(coords[labels,1], d2ys, 'b', lw=3) 
    #     ax.set_ybound(img.shape[0],0)
    #     ax.set_xbound(img.shape[1],0)
In [45]:
show_path_and_der(X, north_bound[:-1])
show_path_and_der(X, north_bound)
c:\Users\louiery\AppData\Local\Continuum\Anaconda\lib\site-packages\numpy\lib\polynomial.py:587: RankWarning: Polyfit may be poorly conditioned
  warnings.warn(msg, RankWarning)
c:\Users\louiery\AppData\Local\Continuum\Anaconda\lib\site-packages\numpy\lib\polynomial.py:587: RankWarning: Polyfit may be poorly conditioned
  warnings.warn(msg, RankWarning)

[  1.06930365e-12  -1.56706172e-09   7.26944248e-07  -2.53057210e-05
  -7.60723893e-02   2.17899060e+01  -1.86671551e+03]
[  7.97533935e-14  -1.44710550e-10   9.30944152e-08  -1.74321336e-05
  -7.36882257e-03   4.41364011e+00  -8.42950169e+02   5.76460830e+04]

In [46]:
show_path_and_der(X, west_bound[:-1])
show_path_and_der(X, west_bound)
c:\Users\louiery\AppData\Local\Continuum\Anaconda\lib\site-packages\numpy\lib\polynomial.py:587: RankWarning: Polyfit may be poorly conditioned
  warnings.warn(msg, RankWarning)
c:\Users\louiery\AppData\Local\Continuum\Anaconda\lib\site-packages\numpy\lib\polynomial.py:587: RankWarning: Polyfit may be poorly conditioned
  warnings.warn(msg, RankWarning)

[  1.06930365e-12  -1.56706172e-09   7.26944248e-07  -2.53057210e-05
  -7.60723893e-02   2.17899060e+01  -1.86671551e+03]
[  4.56575481e-15  -7.34342113e-12   3.99901906e-09  -4.45469719e-07
  -3.82229974e-04   1.64063179e-01  -2.53433005e+01   1.40317464e+03]

In [47]:
show_path_and_der(X, south_bound[:-1])
show_path_and_der(X, south_bound)
c:\Users\louiery\AppData\Local\Continuum\Anaconda\lib\site-packages\numpy\lib\polynomial.py:587: RankWarning: Polyfit may be poorly conditioned
  warnings.warn(msg, RankWarning)
c:\Users\louiery\AppData\Local\Continuum\Anaconda\lib\site-packages\numpy\lib\polynomial.py:587: RankWarning: Polyfit may be poorly conditioned
  warnings.warn(msg, RankWarning)

[  1.06930365e-12  -1.56706172e-09   7.26944248e-07  -2.53057210e-05
  -7.60723893e-02   2.17899060e+01  -1.86671551e+03]
[  1.40141713e-14  -2.46802758e-11   1.53587051e-08  -2.72802762e-06
  -1.18949927e-03   6.79270225e-01  -1.25446024e+02   8.31819923e+03]

In [51]:
def cost_of_path(coords, labels):
    labels = tuple(labels)
    def get_der_coef(coords, labels):
        coefficients = np.polyfit(coords[labels,1], coords[labels,0], len(labels)-1)
        dcoefficients = np.polyder(coefficients)
        return dcoefficients
    initial = get_der_coef(coords, labels[:-1])
    decision = get_der_coef(coords, labels)[1:]
    import scipy
    return scipy.spatial.distance.euclidean(initial, decision)
In [52]:
cost_of_path(X, north_bound)
c:\Users\louiery\AppData\Local\Continuum\Anaconda\lib\site-packages\numpy\lib\polynomial.py:587: RankWarning: Polyfit may be poorly conditioned
  warnings.warn(msg, RankWarning)
c:\Users\louiery\AppData\Local\Continuum\Anaconda\lib\site-packages\numpy\lib\polynomial.py:587: RankWarning: Polyfit may be poorly conditioned
  warnings.warn(msg, RankWarning)

Out[52]:
59519.08081535052
In [53]:
cost_of_path(X, west_bound)
c:\Users\louiery\AppData\Local\Continuum\Anaconda\lib\site-packages\numpy\lib\polynomial.py:587: RankWarning: Polyfit may be poorly conditioned
  warnings.warn(msg, RankWarning)
c:\Users\louiery\AppData\Local\Continuum\Anaconda\lib\site-packages\numpy\lib\polynomial.py:587: RankWarning: Polyfit may be poorly conditioned
  warnings.warn(msg, RankWarning)

Out[53]:
3270.2298416048066
In [54]:
cost_of_path(X, south_bound)
c:\Users\louiery\AppData\Local\Continuum\Anaconda\lib\site-packages\numpy\lib\polynomial.py:587: RankWarning: Polyfit may be poorly conditioned
  warnings.warn(msg, RankWarning)
c:\Users\louiery\AppData\Local\Continuum\Anaconda\lib\site-packages\numpy\lib\polynomial.py:587: RankWarning: Polyfit may be poorly conditioned
  warnings.warn(msg, RankWarning)

Out[54]:
10185.978954457383
In []: